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Abstract 

This  paper  describes  a  model-based  controller  for  the  regulation  of  a  proton  exchange  membrane  (PEM)  fuel  cell.  The  model  accounts 
for  spatial  dependencies  of  voltage,  current,  material  flows,  and  temperatures  in  the  fuel  channel.  Analysis  of  the  process  model  shows  that 
the  effective  gain  of  the  process  undergoes  a  sign  change  in  the  normal  operating  range  of  the  fuel  cell,  indicating  that  it  cannot  be  stabilized 
using  a  linear  controller  with  integral  action.  Consequently,  a  nonlinear  model-predictive-controller  based  on  a  simplified  model  has  been 
developed,  enabling  the  use  of  optimal  control  to  satisfy  power  demands  robustly.  The  models  and  controller  have  been  realized  in  the 
MATLAB  and  SIMULINK  environment.  Initial  results  indicate  improved  performance  and  robustness  when  using  model-based  control  in 
comparison  with  that  obtained  using  an  adaptive  controller. 

©  2004  Elsevier  B.V.  All  rights  reserved. 

Keywords:  PEM  fuel  cell  dynamics;  Model-based  control;  Model  predictive  control;  Transient  model 


1.  Introduction 

Fuel  cells  are  chemical  engines  that  convert  chemical  po¬ 
tential  into  electrical  power.  Since  they  are  not  based  on  tem¬ 
perature  differences,  they  are  not  subjected  to  Carnot’s  limit 
of  efficiency.  In  addition,  common  pollutants  such  as  sulfur 
dioxide  and  nitrous  oxides  are  avoided  since  the  process  does 
not  involve  combustion.  These  advantages,  together  with  the 
reduction  of  greenhouse  gases  and  fuel  consumption  due  to 
higher  efficiencies  and  the  possibility  of  alternative  energy 
sources,  have  generated  enormous  interest  in  fuel  cells  for 
stationary  as  well  as  mobile  applications. 

The  heart  of  the  fuel  cell  system  consists  of  two  elec¬ 
trodes:  the  anode  and  cathode.  The  most  basic  system  uses 
pure  hydrogen  as  fuel,  which  is  oxidized  at  the  anode,  pro¬ 
ducing  electrons  and  protons.  The  electrons  are  released  to 
an  external  circuit,  where  they  can  be  used  to  perform  work, 
while  the  protons  diffuse  through  a  membrane  to  the  cath¬ 
ode.  At  the  cathode,  oxygen  reacts  with  the  electrons  from 
the  external  circuit  and  protons  from  the  anode  reaction, 
forming  water.  Other  systems  are  being  developed,  includ¬ 
ing  the  use  of  methanol  as  fuel  (fed  directly  to  the  fuel  cell 
[10],  or  reformed  externally  to  produce  hydrogen)  and  even 
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using  traditional  fossil  fuels.  Whatever  the  actual  configura¬ 
tion,  the  basic  principles  remain  the  same. 

The  dynamics  of  fuel  cells  are  complex  and  include  the 
mass  transport  of  materials  through  the  membrane  and  to 
and  from  the  electrodes,  reaction  mechanisms  and  rates 
at  the  electrodes,  voltages  and  currents  produced  depend¬ 
ing  on  the  pressures,  temperatures  and  concentrations  at 
the  electrodes,  overpotential  and  Ohmic  losses  of  voltage 
throughout  the  process.  Several  alternative  models  of  dif¬ 
ferent  complexity  have  been  proposed  to  describe  the  per¬ 
formance  of  fuel  cells  under  an  array  of  conditions  (e.g., 
[4,12,8]).  These  models  are  then  used  to  evaluate  optimal 
schemes  of  external  heating,  water  management  [5]  and  fuel 
composition. 

The  dynamic  response  of  fuel  cells  is  important  for  vehic¬ 
ular  applications,  where  power  demands  fluctuate,  and  the 
fuel  cell  does  not  usually  operate  at  the  optimal  steady-state 
designed  by  the  fuel  cell  manufacturer.  A  empirical  model 
for  the  transient  response  of  a  fuel  cell  was  developed  by 
[1].  However,  this  model  is  lumped  and  therefore  does  not 
accurately  portray  the  spatial  dependencies  in  the  system. 
In  a  similar  fashion  Kang  et  al.  [7]  present  an  analysis  of 
a  dynamic  model  for  a  molten-carbonate  fuel  cell  (MCFC), 
where  the  system  is  modeled  as  a  collection  of  first  order 
transfer  functions  with  dead  times.  Pukrushpan  et  al.  [9] 
present  a  nonlinear,  zero-dimensional,  dynamic  model  of  a 
fuel  cell,  but  this  model  is  also  empirical. 
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Nomenclature 

A 

heat  exchange  area  per  unit  length  (cm) 

a 

water  activity 

c 

concentration  at  membrane  surface  (mol/cm3) 

cp 

mass  heat  capacity  ( J/(g  K)) 

d 

channel  height  (cm) 

D 

diffusion  coefficient  in  diffusion  layer  (cm2/s) 

Dq 

intradiffusion  coefficient  of  water  in 
membrane  (cm2/s) 

D* 

diffusion  coefficient  of  water  in 
membrane  (cm2/s) 

e 

membrane  area  per  unit  length  (cm) 

f 

cross-section  of  solid  (cm2) 

F 

Faraday’s  constant  (Col/equivalent) 

h 

channel  width  (cm) 

i 

current  density  (A/cm2) 

A) 

exchange  current  density  (A/cm2) 

k 

heat  conduction  coefficient  (W/(cmK)) 

kc 

condensation  rate  constant  (s_1) 

kp 

water  permeability  (cm2) 

L 

channel  length  (cm) 

M 

molar  flow  (mol/s) 

Win,  dry 

membrane  dry  weight  (g/mol) 

«d 

electro-osmotic  coefficient  of  water  in 
membrane  (molecules/protons) 

P 

pressure  (atm) 

T 

temperature  (K) 

membrane  thickness  (cm) 

u 

convective  heat  transfer  coefficient 
(W/(cm2  K)) 

V 

voltage  (V) 

Greek  symbols 

P 

Water  viscosity  (g/(cm  s)) 

a 

Ratio  of  water  molecules  per  proton 
flux  (molecules/protons) 

s 

Diffusion  layer  thickness  (cm) 

AH 

Enthalpy  of  overall  reaction  (J/mol) 

A^vap 

Enthalpy  of  water  evaporation  (J/mol) 

p 

Density  (g/cm3) 

Pm,  dry 

Dry  membrane  density  (g/cm3) 

G 

Membrane  conductivity  (Ohm-1  cm-1) 

Subscript 

a 

Anode 

avg 

Average 

c 

Cathode 

cool 

Coolant 

g 

Gas 

h2 

Hydrogen 

in 

Inlet 

inf 

Surroundings 

n2 

Nitrogen 

o2 

Oxygen 

oc 

Open  circuit 

s 

Solid 

set 

Set  point 

w 

Water 

Superscript 

1 

Liquid 

sat 

Saturation 

V 

Vapor 

This  paper  focuses  on  the  crucial  issue  of  how  to  control 
the  fuel  cell  to  ensure  acceptable  response  time  for  the  power 
demand,  while  achieving  high  efficiencies  over  the  entire 
operating  range.  To  assist  in  addressing  this  issue.  Section  2 
presents  a  spatial,  time-dependent  model  of  a  fuel  cell.  The 
analysis  of  fuel-cell  controllability  follows,  relying  on  the 
developed  model,  as  well  as  the  synthesis  of  an  adaptive 
controller,  intended  to  account  for  the  observed  sign-change 
in  the  process  static  gain.  Section  3  describes  an  alternative 
approach,  relying  on  a  reduced-order  model,  and  the  use 
of  the  model  in  a  robust  model  predictive  control  (MPC) 
scheme,  which  satisfies  the  power  demand  over  a  wide  range 
of  conditions,  and  is  demonstrated  to  provide  performance 
superior  to  that  of  the  adaptive  scheme. 


2.  Model  formulation  and  solution 

The  model  is  based  on  the  concept  presented  by  [12], 
where  a  fuel  cell  is  modeled  along  its  channel,  as  shown 
in  Fig.  1.  The  model  accounts  for  heat  transfer  between 
the  solid  and  the  two  gas  channels,  and  between  the  solid 
and  cooling  water.  The  water  content  is  modeled,  as  well. 
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Fig.  1.  Schematic  diagram  of  fuel  cell  channel. 
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calculating  the  condensation  and  evaporation,  water  drag 
through  the  membrane,  and  water  generation  at  the  cath¬ 
ode.  The  energy  balance  on  the  solid  is  modeled  dynam¬ 
ically  whereas  all  the  other  equations  are  assumed  to  be 
at  quasi-steady-state  for  a  given  solid  temperature  profile. 
Other  equations  vary  slightly  from  those  presented  by  Yi  and 
Nguyen,  as  well.  Namely,  the  heat  generation  term  is  taken 
from  Berger  [3],  since  this  term  can  be  easily  expressed  in 
terms  of  temperature  (although  this  has  not  been  carried  out 
in  the  present  study). 

The  dynamics  of  the  electrochemistry,  as  well  as  the  tran¬ 
sient  response  of  the  fluids  in  the  channels  (anode  cathode 
and  coolant),  are  assumed  to  be  instantaneous  relative  to  the 
thermal  transient  response  of  the  solid.  Consequently,  con¬ 
servation  equations  that  model  the  entire  system,  with  the 
exception  of  the  fuel  cell  solid  support,  are  formulated  in 
quasi-steady-state  [1],  This  greatly  reduces  the  complexity 
of  the  system  since  it  is  reduced  to  a  one-dimensional  prob- 
lem.For  a  given  local  current  density  (A/cm2),  the  reactants 
are  consumed  as  follows 


d  Mh2(-t) 

dx 

dMp2  (x) 

dx 


—  7(x) 

2  F 

(1) 

h 

—  7(x) 

4  F 

(2) 

The  change  in  the  flow  rate  of  liquid  water  in  each  flow 
channel,  which  is  influenced  solely  by  evaporation  and  con¬ 
densation,  is  given  as 


kchd 

dx  Rl\.(x) 

k  —  a,c 


£fkMa(x) 


Pk  ~  Pf(Tk) 


(3) 


This  implies  that  the  rate  of  condensation/evaporation  is 
proportional  to  the  difference  between  the  partial  pressure 
of  water  and  the  saturation  pressure  at  the  local  tempera¬ 
ture.  Note  that  in  the  event  of  no  liquid  water  and  a  partial 
pressure  lower  than  the  saturation  pressure,  this  equation 
is  not  valid.  (See  Appendix  C  on  how  this  is  dealt  with 
during  integration).  The  water  vapor  in  the  flow  channels 
is  affected  by  a  number  of  mechanisms:  (a)  water  vapor 
is  generated  at  the  cathode  by  the  reaction  of  oxygen  with 
the  proton  and  electron;  (b)  water  vapor  diffuses  through 
the  membrane;  (c)  water  vapor  is  dragged  through  the 
membrane  by  migrating  protons  (from  the  anode  to  the 
cathode);  (d)  liquid  water  can  condense  and  evaporate  based 
on  the  partial  pressure  of  water  and  the  saturated  pressure 
(temperature-dependent).  These  mechanisms  are  included 
in  the  following  equation,  which  models  the  anode  flow 
channel 


dMl  (x)  d M]  (x)  ha(x) 

- = - 7(x) 

dx  dx  F 


(4) 


where  the  first  term  on  the  right  is  the  condensation  or 
evaporation  of  the  water  and  the  second  is  the  net  migration 
of  water  across  the  membrane.  This  net  migration  is  the 


combined  effect  of  diffusion  by  concentration  and  pressure 
gradients  along  the  membrane  and  water  molecules  dragged 
across  the  membrane  by  the  current.  The  ratio  of  water 
molecules  per  proton,  a,  is  given  by 


kp  F  d  I\ 


F  n*dcw 

a  =  n  d - D - cw 

/(x)  dy  ii  7(x)  dy 

F  ^w,c  £*w,a 

ss  rid - D  - 

7(x)  tm 

(Cw,a  "T  C\v,c)  kp  F  Py/  c  Pv/,a 
2  fi  7(x)  tm 


(5) 


The  water  vapor  content  of  the  cathode  is  modeled  using 
the  following  equation,  which  includes  the  generation  of 
water  by  reaction 


dM;,c(x) 

dx 


d  M, 


w,c(a)  h  ha(x) 

~a - + 

dx  2  F  F 


(6) 


The  local  temperatures  in  the  anode,  cathode  and  coolant 
channels  are  affected  by  heat  transfer  between  the  mass  sur¬ 
face  and  the  fluid 


d Tk (x)  =  UgAg{Ts(x)  -  Tk(x)} 
dx  Y.,  CpjMj(x) 


drcooi(x)  _  f/cooiAcooi{rs(x)  -  rcooi(x)} 

dx  Cpw  A7coo| 


(8) 


The  cell  voltage  is  inversely  proportional  to  the  current  den¬ 
sity,  and  related  to  the  partial  pressures  of  the  species,  using 
the  Nernst  and  Tafel  equations  [2] 


keel  I  —  k  >e  + 


RT 

IF 


In 


/[%bW-(i|) 


POj.bix)  ~ 

(  SI(x)  \ 

\4FDo2) 

0.5 

^20,bW  +  ( 

SI(x)  ) 

-FDn2o  J 

\ 


V 


/ 


RT]n  I  7(x)  \  _  7(x)fm 

F  °'"w 


This  equation  is  derived  as  shown  in  Appendix  A.  Thus, 
the  effective  cell  voltage  is  the  thermodynamically  reversible 
cell  voltage  minus  the  Ohmic,  activation  and  concentration 
overpotentials.  The  partial  pressures  are  simply  the  molar 
ratios  multiplied  by  the  electrode  pressure 

Mi 

pi  =  ^frP  dO) 


The  following  are  empirical  expressions  of  the  various 
constants  used  in  the  model  equations  (following  [11]) 


Oin  =  (o.00514^^cw,fl 
V  Pdry 


0.00326) 


exp 
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Ow.k- 


w, 


Pdry  (0.043  +  17.8a*  -  39.85^  +  36a^) 


m,dry 


ajc  <  1  k  =  a,  c 
5^(14 +1.4(a*-l)) 

™m,dry 

a*  >  1  k  =  a,  c 


(12) 


M. 


ak  = 


w.k 


pk 


E/  CpjMiix)  P sat (7*  -  273) 


,  k  =  a,  c 


(13) 


p..t(7)  =  j  q—2.  18+2.95  x  10_2r— 9.18xl0_5r2+1.44x  10_7T3 

(14) 

__  |  0.0049  +  2.024aa  -  4.53a2  +  4.09a3  «a  <  1 
”d  _  {  1.59  +  0. 159(aa  -  1)  aa  >  1 

(15) 

D*=ndD„  exp  (2416  (+-!))  (16) 


A//vap(7)=45,070-41.947’  +  3.44  x  10“3P2 

+  2.54810“6r3  -  8.98  x  10“107’4  (17) 


The  heart  of  the  model  is  a  spatial,  time -dependent  partial 
differential  equation  describing  the  local  temperature  of  the 
solid  cell  support 


dTs  d  2Ts  U„A„ 

psCPs =ks—{  +  +  TC-  2 Ts) 

dt  oxz  j 

(7cool2lcool  61  ( AH  , 

+  (Tcool  -  rs)  -  -  [  —  +  Ken  )  m 


+  —  AHmp(T) 


d Mi  Ax)  dM'(x) 


d* 


+ 


djc 


(18) 


with  boundary  conditions 
=  UC(TS  -  7-nf) 


,  97s 
k  — 
dx 


dTs 
k  — 
dx 


x=0 


=  —UC(TS  -  7-nf) 


x=L 


(19) 


Evidently,  the  channel  solid  temperature  is  affected  by:  (1) 
heat  transfer  by  conduction,  (2)  heat  transfer  by  convection 
with  the  flow  and  coolant  channels,  (3)  heat  generation  by 
the  reaction  (total  enthalpy  change  of  the  reaction  (water 
formation)  minus  the  external  power  [3])  and  (4)  heat  of 
condensation/evaporation  of  water  in  the  flow  channels.  The 
boundary  conditions,  shown  in  Eq.  (19),  reflect  heat  losses 
to  the  surroundings  from  the  edges.  It  is  assumed  that  there 
is  no  heat  loss  along  the  channel,  since  the  channel  is  part 
of  a  large,  symmetrical  proton  exchange  membrane  (PEM) 
stack. 


2.1.  Solution  procedure 

Eq.  (18)  is  integrated  using  the  Crank-Nicholson  method, 
in  which  the  partial  differential  equation  is  discretized,  re¬ 
sulting  in  a  system  of  algebraic  equations  that  defines  the 
dependence  between  the  solid  temperature  profiles  for  con¬ 
secutive  times.  For  each  time-step  in  the  integration,  the 
quasi-steady  state  ODEs  (Eqs.  ( 1)— (8))  are  integrated,  and 
the  surface  temperature  profile  is  calculated  for  the  next 
time-step.  Full  details  of  the  Crank-Nicholson  method  and 
calculations  are  presented  in  Appendix  D. 

During  the  integration  of  the  ODEs  the  cell  voltage  is 
set.  If  the  system  is  modeled  for  constant  voltage,  the  aver¬ 
age  current  density  is  calculated  for  every  time  step  of  the 
Crank-Nicholson  integration  by  integrating  the  values  of  I(x) 
over  x  and  dividing  by  the  length  of  the  channel.  However, 
if  the  system  has  a  set  current  density,  the  voltage  is  guessed 
and  iteratively  manipulated  (using  the  secant  method)  un¬ 
til  the  calculated  current  density  converges  to  the  set  value. 
This  procedure  is  problematic,  since  there  is  no  guarantee 
that  the  system  (with  specific  flow  rates,  concentrations  and 
temperatures)  can  attain  the  desired  current  density.  A  better 
alternative,  would  be  to  specify  an  external  load  resistance, 
Rext ,  and  then  match  the  estimated  Vceii  with  the  calculated 
value,  Vcalc  =  7aVg,caic  x  Rext-  This  equation  can  always  be 
fulfilled  since  it  is  based  on  physical  behaviors,  as  opposed 
to  arbitrary  guesses.  This  has  not  been  implemented  in  this 
paper,  to  allow  for  comparison  with  results  in  the  literature. 

The  data  used  for  the  simulation  is  based  on  the  study  by 
Trung  and  Nguyen  (1998),  and  is  presented  in  Table  1. 

The  model  formulated  above  constitutes  a  set  of  implicit 
equations,  since,  as  shown  in  Appendix  D,  the  equations 
depend  nonlinearly  on  the  solid  temperature.  To  reduce  the 
computation  time  inherent  in  the  implicit  formulation,  a 
partially-explicit  approximation  was  developed  as  an  alter¬ 
native,  where  the  nonlinear  expressions  of  the  solid  temper¬ 
ature  are  used  explicitly,  leading  to  a  linear  system  for  each 
time  step.  To  compare  the  difference  between  the  solutions 
obtained  using  the  original  implicit  formulation  with  a  par¬ 
tially  explicit  approximation,  the  system  transient  is  com¬ 
puted  starting  from  a  uniform  solid  temperature  of  343  K 
using  the  data  from  the  base  case.  The  evolution  of  the  av¬ 
erage  temperature  and  power  density  with  time  for  varying 
time-step  sizes  are  presented  for  both  methods  in  Fig.  2,  not¬ 
ing  that  the  system  settles  well  inside  20  s.  As  can  be  seen, 
for  a  step  size  of  0.1  s,  there  is  no  difference  in  the  results 
obtained  from  the  two  methods.  At  0.5  s  there  is  a  small  dif¬ 
ference  between  the  two,  with  the  implicit  method  repeating 
the  results  for  time  step  of  0.1  s.  Since  the  implicit  method 
requires  the  solution  of  nonlinear  equations,  it  is  more  com¬ 
putationally  demanding,  and,  for  small  step  sizes,  unneces¬ 
sary. 

The  spatial  dependence  of  the  variables  at  the  end  of  the 
simulation  time  (i.e.  at  steady  state)  can  be  seen  in  Fig.  3. 
The  plot  shows  the  computed  values  indicated  by  crosses, 
with  a  cubic  spline  used  to  indicate  the  continuous  profile. 
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Fig.  2.  Simulation  of  system  starting  from  uniform  temperature  of  70  °C.  The  evolution  of  the  average  temperature  and  average  power  density  with  time 
using  explicit  and  implicit  Crank-Nicholson  for  different  time  steps. 


(a)  (b) 


o 


Fig.  3.  Spatial  dependencies  of  variables  in  channel  direction:  (a)  solid  temperature;  (b)  current  density;  (c)  cathode  water  vapor  flow;  (d)  cathode  water 
liquid  flow. 
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Table  1 

Data  for  base  case  presented  by  Trung  and  Nguyen 


Variable 

Value 

Physical  data 

Acod  (cm) 

0.4 

Ag  (cm) 

0.4 

d  (cm) 

0.1 

ae  (cm) 

0.1 

’/(cm2) 

0.0106 

h  (cm) 

0.1 

L  (cm) 

10 

tm  (cm) 

0.01275 

aDw  (crrr/s) 

10”6 

“.Do,  (cm2/s) 

10~6 

aDH,  (crrr/s) 

10~6 

3  8  (cm) 

0.001 

Pm, dry  (g/cm3) 

2 

aps  (g/cm3) 

2 

Physical  parameters 

/o  (A/cm2) 

0.01 

Vac  (V) 

1.1 

Do  (cm2/s) 

5.5  x  10~7 

h  (s_l) 

100 

kv  (cm2) 

1.58  x  10^14 

p  (g/cms) 

3.565  x  10“3 

ks  (W/cm2K) 

0.005 

t/c  (W/cm2K) 

0.025 

t/g  (W/cm2K) 

0.025 

£/w  (W/cm2K) 

0.025 

Wm,dry 

1100 

AH  (J/mol) 

-2.41  x  1()5 

aCps  (J/gK) 

1 

Operating  values  for  the  fuel  cell 

MH,.  in  (mol/s) 

1.14  x  10"5 

Mv  ■ 

Saturated 

Ml,a.in  (mol/s) 

0 

T,,in  (K) 

353 

M0,. in  (mol/s) 

5.7  x  10-6 

K.a, in  (mol/s) 

0 

K,a, in  (mol/s) 

0 

7c, in  (K) 

353 

Pa  (atm) 

1 

Pc  (atm) 

1 

Tint  (K) 

343 

/avg  (A/cm2) 

1.1 

a  Denotes  data  that  were  estimated  for  this  paper. 


Clearly,  the  distributed  variables  feature  significant  special 
variation  along  the  channel  length. 

3.  Controllability  analysis 

To  obtain  smaller,  lighter  and  cheaper  fuel  cells,  they 
should  be  designed  to  operate  at  the  highest  possible  power 
density.  In  this  regard,  Fig.  4  displays  the  power-current 
relationship  for  the  fuel  cell  model,  operated  at  constant 
flowrates  of  the  feed  streams.  A  similar  result  is  obtained 
when  using  a  constant  fuel  utilization  factor. 

To  provide  guidance  for  the  analysis  of  the  system  dy¬ 
namics  and  its  controllability,  a  series  of  step  tests  were  per- 


Fig.  4.  Power-current  curve  for  base  case. 


formed  on  the  base  case  system.  These  step  tests  are  helpful 
in  determining  the  influence  of  various  variables  on  the  con¬ 
trolled  output.  In  particular,  it  is  of  interest  to  examine  the 
influence  of  the  cell  voltage  on  the  power  density,  as  shown 
in  Fig.  5.  Note  that  the  response  of  the  power  density  to 
changes  in  the  cell  voltage  consists  of  a  large  initial  jump, 
followed  by  a  gradual  change.  The  initial  response  is  the 
electrochemical  dependence,  which  is  instantaneous,  having 
assumed  quasi  steady  state.  The  gradual  change  is  due  to 
the  slow  transient  response  of  the  solid  temperature  profile, 
which  affects  the  electrochemical  behavior  along  with  the 
quasi-steady-state  conditions  in  the  channels.  Furthermore, 
as  can  be  seen  in  Fig.  5,  the  response  of  the  power  to  positive 
step  changes  in  the  cell  voltage  depends  on  the  operating 
level  selected,  with  the  net  change  of  the  power  density  be¬ 
ing  either  positive  at  low  voltages  or  negative  at  high  ones. 

These  changes  in  the  static  gain  ( AP/AV)  are  reflected 
in  Fig.  6  where  the  static  gain  is  plotted  as  a  function  of 
average  current  density  for  two  sets  of  conditions:  the  base 
case  and  slightly  modified  conditions:  namely,  reactant  and 
oxidant  inlet  temperatures  of  70  and  60  °C,  respectively, 
and  a  hydrogen  and  oxygen  inlet  flowrates  of  1.71  x  10-5 
and  8.55  x  10-6  mol/s,  respectively.  As  can  be  seen  in 
Fig.  6,  the  static  gain  changes  sign  for  the  base  case  at  a 
current  density  of  about  0.9  A/cm2.  On  the  other  hand,  for 
the  modified  conditions,  the  sign  change  occurs  at  a  higher 
current  density-around  1.3  A/cm2.  Alternatively,  if  the  cur¬ 
rent  density  is  used  as  the  control  variable,  the  gain  be¬ 
tween  the  power  density  and  the  average  current  density 
(AP/A/avg)  must  be  evaluated,  as  shown  in  Fig.  7.  Once 
again,  a  sign  change  is  evident.  For  the  base  case,  this  occurs 
at  roughly  1.1  A/cm2.  Again,  for  the  modified  conditions 
the  sign  change  occurs  at  a  different  operating  point-about 
0.85  A/cm2.  Note  that  the  static  gain  changes  both  in  mag¬ 
nitude,  but  more  importantly,  in  sign  as  well.  In  addition, 
other  variables  change  the  location  of  the  sign  change.  This 
creates  severe  problems  for  a  control  system,  since  the  di- 
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Fig.  5.  Step  tests  for  base  case.  In  each  case,  a  step  change  of  0.05  V  is  implemented,  from  operating  levels:  (a)  Vcen  =  0.15  V,  (b)  Vceii  =  0.40  V,  (c) 
Vceii  =  0.60  V  and  (d)  Vceii  =  0.85  V. 


rection  of  the  power’s  change  in  response  to  the  control  vari¬ 
able  (voltage  or  current  density)  is  not  constant  and,  over 
a  large  part  of  the  operating  space,  difficult  to  predict  ac¬ 
curately.  The  sign-change  in  static  gain  precludes  the  use 


Fig.  6.  Static  gain  of  power  density  to  changes  in  cell  voltage  for  base 
case  and  modified  conditions. 


of  a  fixed-gain  controller  with  integral  action,  since  such  a 
controller  cannot  be  stabilized.  Two  alternative  options  that 
can  be  made  to  work  are  adaptive  control  and  nonlinear 
model-based  control. 


Fig.  7.  Static  gain  of  power  density  to  changes  in  average  current  density 
for  base  case. 
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A  simple  adaptive  controller  attempts  to  compensate  for 
the  changing  process  gain  by  adapting  itself  on  the  fly.  This 
is  difficult  to  design  since  the  actual  behavior  of  the  process 
is  dependent  on  a  number  of  variables  (flowrates,  compo¬ 
sitions,  temperatures,  etc.)  as  demonstrated  in  Fig.  6.  Thus, 
an  adaptive  controller  that  does  not  account  for  all  of  these 
inputs  will  experience  difficulties,  and  consequently,  will 
have  limited  performance  and  robustness.  Alternatively,  a 
fixed-gain  controller  can  be  implemented  to  control  the  fuel 
cell  in  ranges  of  operating  conditions  where  a  sign  change 
in  the  static  gain  is  not  expected.  There  are  two  problems 
with  this  approach: 

1.  This  kind  of  system  is  constantly  in  danger  of  straying 
into  the  “instability  zone,”  and,  once  the  system  gets 
there,  a  controller  including  integral  action  and  previ¬ 
ously  stable  will  become  unstable  when  the  process  gain 
changes  direction. 

2.  The  operating  conditions  chosen  to  ensure  stable  opera¬ 
tion  are  unlikely  to  be  optimal  in  terms  of  fuel  efficiency 
and  power  density. 

A  control  system  that  relies  on  a  nonlinear  model  can 
avoid  both  problems.  The  model  will  indicate  the  actual 
influence  the  control  variable  will  have  on  the  system. 
This  is  the  motivation  to  formulate  a  model  that  predicts 
the  behavior  with  sufficient  accuracy  while  being  simple 
enough  to  enable  its  use  in  an  on-line  optimal  control 
scheme. 

3.1.  Adaptive  control 

When  designing  a  control  system,  a  number  of  re¬ 
quirements  need  to  be  satisfied.  The  control  system  must 
eliminate  disturbances  and  track  changes  in  the  set  point 
accurately  and  with  sufficient  speed.  This  must  be  done  for 
as  wide  a  set  of  conditions  as  possible.  Also,  it  is  desirable 
that  the  control  system  be  as  simple  as  possible.  Therefore, 
the  design  of  a  simple  adaptive  controller  was  attempted 


first,  with  the  experience  serving  also  as  an  illustration  of 
some  of  the  phenomena  discussed  earlier. 

Examining  Fig.  7  reveals  the  near  linear  dependency  of 
the  static  gain  on  the  current  density.  Therefore,  for  control 
purposes,  a  standard  PI  controller  could  use  an  adaptive  gain 
inversely  proportional  to  the  variable  process  static  gain 


G(/aVg)  1.5s(j+  1)  (20) 

G(/avg)  =  — 1.333/avg  +  1.466 

To  smooth  out  the  otherwise  potentially  abrupt  changes 
in  the  controller  gain,  the  gain  value  is  adjusted  using  a 
low  pass  filter  with  a  time  constant  of  1  s.  Fig.  8  presents 
the  performance  obtained  with  this  adaptive  controller.  Note 
that  it  performs  nicely  for  the  base  case  (for  which  it  was 
designed),  converging  within  roughly  6  s  for  positive  and 
negative  step  changes  in  the  power  density  set  point.  How¬ 
ever,  for  the  modified  conditions,  it  cannot  reach  the  set 
point  of  0.55  W/cm2  since  the  controller  changes  the  gain’s 
sign  at  1.1  A/cm2,  instead  of  at  0.85  A/cm2.  Note  that  the 
response  for  the  decrease  is  more  aggressive,  with  over¬ 
shoot  because  of  the  inaccurate  value  of  the  static  gain  that 
the  controller  uses.  Full  design  of  a  standard  controller  re¬ 
quires  that  the  uncertainty  of  the  actual  process  gain  be 
addressed. 


4.  Model  predictive  control 

As  has  been  shown,  the  fuel  cell  presents  a  number  of  con¬ 
trol  problems,  the  most  significant  of  which  is  the  nonlinear¬ 
ity  in  the  vicinity  of  the  peak  power  density.  Furthermore,  it 
is  of  interest  to  ensure  efficient  operation  during  transients 
and  when  extreme  load  changes  are  imposed.  Consequently, 
it  is  proposed  to  tackle  these  problems  using  a  nonlinear 
model  predictive  control  framework. 

Model  predictive  control  is  part  of  a  family  of  optimiza¬ 
tion-based  control  methods,  which  are  based  on  on-line 


Fig.  8.  Control  of  set  point  step  changes  for  base  case  (left)  and  modified  conditions  (right). 
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optimization  of  future  control  moves.  Using  a  process  model, 
the  optimizer  predicts  the  effect  of  past  inputs  on  future  out¬ 
puts.  Then,  using  the  same  model,  it  computes  a  sequence 
of  future  control  moves,  such  that  an  objective  function,  in¬ 
cluding  penalties  on  the  trajectory  of  predicted  tracking  er¬ 
ror,  is  minimized.  The  first  of  the  future  control  moves  is 
implemented,  and  the  entire  optimization  is  repeated  from 
the  next  step  on,  and  so  on,  ad  infinitum.  Feedback  is  used 
to  compensate  for  the  model’s  inaccuracies  and  to  ensure 
convergence.  For  a  detailed  account  of  nonlinear  MPC,  the 
reader  is  referred  to  the  excellent  review  by  Henson  [6], 

The  definition  of  the  optimization  problem  offers  enor¬ 
mous  advantages  for  advanced  control.  Typically,  the  prob¬ 
lem  formulation  includes  the  satisfaction  of  the  set  points 
for  the  different  output  variables  (least  squared  errors).  The 
optimization  problem  uses  a  model  of  the  system  to  pre¬ 
dict  the  output  variables  values  and  compares  them  with 
the  set  points.  Hence,  the  performance  of  the  controller  is 
influenced  by  the  accuracy  of  the  model.  In  addition,  the 
optimization  problem  commonly  includes  different  degrees 
of  penalty  for  changes  in  the  control  variables.  This  en¬ 
courages  the  controller  to  avoid  excessive  changes  of  the 
most  “expensive”  control  variables.  The  advantage  of  the 
flexibility  of  the  optimization  problem  is  apparent.  Differ¬ 
ent  degrees  of  importance  can  be  attributed  to  set  point 
satisfaction,  control  variable  step  size  and  absolute  values 
(influencing  stability  and  controller  behavior),  and  other 
values  of  varying  importance.  For  instance,  the  model’s 
prediction  of  overall  efficiency  can  be  taken  into  considera¬ 
tion.  The  optimization  can  be  unconstrained  or  may  include 
hard  constraints  (for  example  safety  limits  or  downstream 
demands  on  the  process  being  controlled). 

The  first  step  in  designing  an  MPC  system  is  the  derivation 
of  a  model  that  the  controller  will  use  for  the  optimization. 
This  model  should  be  as  accurate  as  possible,  while  being 
simple  enough  to  allow  for  repeated  calculations  during  the 
optimization. 


4.1.  Reduced  model  -  multiple  CSTRs 


This  simplified  model  neglects  the  spatial  dependence  of 
the  variables  in  the  flow  channels.  In  essence,  this  assumes 
CSTR  conditions  of  the  reactants  and  oxidant,  and  of  the 
coolant,  as  well.  Therefore,  differential  Eqs.  ( 1)— (8)  reduce 
to  the  following  algebraic  equations 
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Similarly,  the  partial  differential  Eq.  (18)  reduces  to  the 
ordinary  differential  equation 

d  Ts  Uga  Ucooib 

PsCps-f  =  ^-(T.d  +  Tc-  2 Ts)  +  (Tcooi  -  Ts) 


d  t  f 
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The  current  density  is  the  solution  of  Eq.  (9).  The  solution 
method  is  to  solve  all  the  algebraic  equations  for  each  time 
step  of  the  integration  of  Eq.  (28).  The  solution  of  a  and 
the  current  density  is  cascaded  in  the  following  fashion,  as 
shown  in  Fig.  9: 

1.  Receive  the  cell  voltage  Ucen. 

2.  Guess  the  cell  current  density  I. 

3.  Calculate  the  molar  flows  of  oxygen  and  hydrogen. 

4.  Guess  ct\ . 

5.  Use  the  solid  temperature  to  calculate  the  saturated  pres¬ 
sure  in  the  anode  and  cathode  (this  simplifies  the  solu¬ 
tion). 

6.  Calculate  the  molar  flows  of  the  water  vapors  (see  be¬ 
low). 

7.  Calculate  the  temperatures  of  the  flow  channels. 

8.  Calculate  a  and  converge  (step  4). 

9.  Calculate  the  voltage  and  compare  to  set  voltage 

10.  Change  current  accordingly  (step  2). 

The  calculation  of  the  water  content  in  the  anode  and  cath¬ 
ode  for  a  given  current  density  is  somewhat  tricky.  Eqs.  (23) 
and  (24)  are  appropriate  if  the  liquid  and  vapor  in  the  chan¬ 
nel  are  in  equilibrium.  This  means  that  liquid  water  remains 
in  the  channel,  either  by  condensation  or  by  only  partial 
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Fig.  9.  Flowchart  for  the  solution  of  the  reduced  model. 


evaporation  of  feed  liquid  water.  Rearranging  Eq.  (24)  re¬ 
sults  in  the  following  expression 
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This  is  easily  solved,  giving  two  solutions  for  the  water 
vapor  content.  For  the  solution  to  be  physically  meaningful, 
the  solution  must  be  greater  than  zero  and  noncomplex.  The 
liquid  water  is  then  calculated  using  Eq.  (23).  If  there  is  no 
physical  solution  for  Eq.  (29)  it  means  that  liquid  water  is 
not  in  equilibrium  with  the  vapor.  In  other  words,  the  liquid 
flow  out  of  the  channel  is  zero  and  the  water  vapor  is  the  sum 
of  the  inlet  flow  water  vapor  and  liquid  water  subtracting 
the  water  migrating 
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If  the  result  of  this  calculation  is  negative  vapor  flow,  it 
means  that  the  guessed  values  of  a  and  I  are  not  feasible,  so 
for  the  sake  of  the  calculation,  is  taken  as  zero.  The 


migrating  water  is  simply  the  difference  between  the  inlet 
flows  of  water  (liquid  and  vapor)  and  the  outlet  flows 

migrate  =  +  M^a  in  -  Mvw  a  -  M]w  a  (31) 

This  value  is  then  used  in  the  calculation  for  the  water 
content  of  the  cathode.  In  a  similar  fashion  to  the  anode, 
Eq.  (26)  is  rearranged  to  give 
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whose  solution  is  used  in  Eq.  (23)  to  calculate  the  liquid 
water  flow  in  the  cathode.  Once  again,  if  neither  result  of 
Eq.  (32)  is  physically  feasible,  the  liquid  flow  is  zero,  and 
the  vapor  flow  is  calculated  using 

,  hL 

Ml,a  =  K,a,in  +  Ml,a,in  +  “is™16  +  ^ (33) 

<a  =  0 

Note  that  this  includes  water  vapor  generated  by  the  re¬ 
action. 

At  high  current  densities,  there  are  problems  converging 
the  current,  since  it  assumed  that  there  is  no  pressure  drop 
in  the  channel.  Obviously,  if  gas  molecules  are  depleted 
(by  reaction  and  condensation),  the  pressure  will  drop,  as 
well.  This  causes  problems  since  if  the  total  pressure  drops, 
the  partial  pressures  drop,  as  well.  Otherwise,  the  partial 
pressures  are  dependent  solely  on  the  ratios  between  the 
components.  Therefore,  the  partial  pressures  are  calculated 
using 

Mi 

Pi  =  ot  (34) 

£M i 


Assuming  ideal  gases,  we  can  use 
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to  relate  the  pressure  to  the  temperature  and  mole  flows  in 
the  channel.  This  results  in 
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(36) 


which  expresses  the  dependence  of  the  partial  pressures  of 
the  components  on  the  mole  flows  and  the  temperature  in 
the  channels.  The  reference  values  are  those  at  the  channel 
inlet. 

This  model  assumes  uniform  conditions  at  all  points  in  the 
channel.  If  this  leads  to  unsatisfactory  accuracy,  dividing  the 
system  into  a  number  of  consecutive  CSTRs  can  incorporate 
a  certain  amount  of  spatial  dependence.  In  this  case,  the 
output  of  one  is  the  input  of  the  next.  The  assumptions  used 
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Fig.  10.  Comparison  of  the  full  model  between  the  multiple  CSTR  model  using  100  CSTRs  with  no  heat  conduction  -  the  average  current  density  vs. 
time  and  the  current  density  profile  at  steady  state. 


in  this  model  are: 

1.  There  is  no  heat  transfer  by  conduction  between  two 
neighboring  cells. 

2.  Heat  is  only  lost  to  the  surroundings  from  the  two  outer¬ 
most  CSTRs. 

3.  The  inlet  temperatures  are  used  to  calculate  the  saturated 
pressures  and  the  physical  parameters  in  the  channels. 
This  is  done  because  the  water  flow  depends  on  the  chan¬ 
nel  temperatures  (condensation).  The  current,  in  turn,  de¬ 
pends  on  the  water  content,  and  the  temperatures  depend 
on  the  current.  If  the  inlet  temperatures  weren’t  used,  the 
result  would  be  a  nonlinear  equation. 

4.1.1.  Results 

The  following  results  display  the  reduced  model’s  ability 
to  predict  the  behavior  of  the  full-order  model.  First,  it  is 
necessary  to  ensure  that  the  CSTR  model  is  appropriate  by 
comparing  the  full  model  with  a  large  number  of  CSTRs.  We 
expect  that  the  two  models  predict  identical  profiles  of  the 
temperatures  and  concentrations  for  any  given  solid  temper¬ 


Fig.  1 1 .  Open  loop  response  to  a  step  change  in  the  current  density  from 
0.5  to  0.8  A/cm2. 


ature  profile.  When  considering  the  transient  behavior  over 
time,  it  is  important  to  realize  that  the  multiple  CSTR  model 
will  not  agree  entirely  with  the  full  model,  since  it  does  not 
account  for  heat  exchange  by  conduction  like  the  full  model. 
The  asymptotic  case  is  if  the  heat  transfer  by  conduction  co¬ 
efficient  is  zero  or  negligible.  The  results  of  the  full  model 
and  the  multiple  CSTR  are  presented  in  Fig.  10.  As  can  be 
seen,  the  two  models  are  in  excellent  agreement. 

It  is  important  for  the  model  to  accurately  predict  the  tran¬ 
sient,  open  loop  behavior  of  a  fuel  cell  under  varying  condi¬ 
tions,  since  the  performance  of  a  control  system  relying  on 
this  model  depends  directly  on  the  model’s  fidelity.  Fig.  11 
compares  open  loop  response  of  the  system  to  a  step  change 
in  the  current  density  from  an  initial  value  of  0.5  A/cm2  to 
0.8  A/cm2  using  the  full  model  and  the  reduced  model  using 
varying  numbers  of  CSTRs.  Fig.  12  shows  the  same  behav¬ 
ior  for  a  current  density  step  change  from  1.1  to  1.4  A/cm2. 

As  can  be  seen,  the  prediction  capabilities  of  the  reduced 
model  are  greater  at  low  current  densities  than  at  high  densi¬ 
ties.  This  makes  sense  since  the  spatial  dependencies  exhibit 
significantly  higher  variation  at  higher  values  of  current 


Fig.  12.  Open  loop  response  to  a  step  change  in  the  current  density  from 
1.1  to  1.4  A/cm2. 
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density,  and  the  reduced  model  neglects  this  variation.  In 
Fig.  12,  one  can  see  that  the  10-CSTR  model,  while  failing 
to  predict  the  value  accurately,  still  retains  the  characteristic 
behavior  nicely,  whereas  the  single  CSTR  model  accurately 
predicts  the  steady  state  change  in  power  density.  This  is 
important  to  remember  when  discussing  model-based  con¬ 
trol,  since  the  quality  of  the  controller  depends  directly  on 
the  quality  of  the  model. 

4.2.  Amphlett-based  reduced  model 

The  solution  of  the  previously  mentioned  model  for  a  20  s 
horizon  occupies  2  s  of  CPU  time  on  a  2.4  GHz  computer. 
This  is  clearly  too  slow  for  practical  purposes,  since  this 
calculation  is  at  the  heart  of  the  MPC  objective  function. 
Moreover,  there  is  no  way  to  calculate  an  analytical  gradient 
for  the  optimizer  forcing  the  optimizer  to  calculate  one  nu¬ 
merically,  at  each  iteration.  Therefore,  a  new  method  is  in¬ 
troduced,  based  on  Amphlett  [1],  where  the  current  density 
is  set,  and  the  cell  voltage  is  calculated.  Amphlett  used  ma¬ 
terial  and  energy  balances  coupled  with  a  parametric  model, 
which  forecasts  the  electrochemical  behavior  to  model  the 
dynamics  of  a  fuel  cell.  This  paper  uses  Eqs.  (9)— (16)  to 
predict  the  electrochemical  conditions.  As  the  model  sup¬ 
ports  no  spatial  variation,  it  is  good  only  for  a  single  CSTR 
description  of  a  fuel  cell. 

In  this  method,  the  current  density  is  set  and  the  only  it¬ 
erations  performed  provide  the  solution  of  the  water  content 
as  demonstrated  previously.  In  this  formulation,  the  input 
variable  is  the  current  density,  and  once  the  water  content  is 
established  for  that  current  density,  the  voltage  can  be  ex¬ 
plicitly  calculated  using  Eq.  (9). 

4.2.1.  Results 

Fig.  13  demonstrates  that  the  Amphlett-based  model  is 
identical  to  the  1-cell  CSTR  model.  In  this  example  the 
current  density  is  changed  from  0.9  to  1.1  A/cm2. 

As  can  be  seen,  the  two  methods  are  in  complete  agree¬ 
ment.  This  is  hardly  surprising  since  the  two  methods  rely  on 


identical  equations  and  assumptions,  differing  only  in  nume¬ 
rical  approach.  Note  that  Amphlett’ s  approach  is  significa¬ 
ntly  faster,  but  cannot  be  expanded  to  more  than  one  CSTR. 

4.3.  The  optimization  problem 


The  objective  function  for  the  MPC  controller  is  the  min¬ 
imization  of  the  sum  of  squared  errors  between  the  desired 
set  point  and  the  actual  trajectory  of  the  power  output,  with 
an  additional  penalty  imposed  on  rapid  changes  in  the  ma¬ 
nipulated  variables 


f(u)  = 


f[ 


w(t)(P(u,  t)  -  pset(t)Y 


d  t 


(37) 


The  weight  functions  W  and  .S',  are  used  to  increase  the 
importance  of  specific  variables  at  given  instances.  For  ex¬ 
ample,  the  weights  may  increase  over  time  to  ensure  rapid 
convergence  with  no  offset.  The  function  is  discretized  over 
time,  obtaining  the  following  algebraic  function 


f(u)  =  ( P(u )  -  Pset)T  W{P(u)  -  Pset)  +  y^  dujT S,dUj 

i 

(38) 


The  vector  P(u)  is  the  actual  value  of  the  power  density 
at  the  different  time  steps  in  the  prediction  horizon,  while 
element  i  of  vector  d»  is  the  value  of  u  at  time  step  i  minus 
its  value  at  time  step  ;  —  1 . 

Note  that  the  actual  variables  in  the  optimization  are  the 
changes  in  values  from  each  time  step  to  the  next.  This 
means  that  the  value  of  u  at  time  step  i  is  simply  the  initial 
value  of  u  plus  all  the  values  of  du  up  to  time  step  i. 


4.3.1.  Results 

Fig.  14  shows  the  performance  obtained  with  MPC  for  the 
modified  conditions  of  the  base  case  as  in  Fig.  8  (using  adap- 


Fig.  13.  Comparison  of  Amphlett  based  model  and  1  CSTR. 


Fig.  14.  MPC  response  for  modified  conditions  of  base  case. 
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Fig.  15.  MPC  response  for  modified  conditions  of  base  case  to  power 
density  of  0.54  W/crrr. 

five  control).  As  before,  the  current  density  cannot  satisfy  the 
increase  in  power  demand,  but  this  time,  the  controller  does 
not  experience  loss  of  stability.  Furthermore,  the  response 
of  MPC  in  tracking  a  demand  to  decrease  the  power  output 
performs  significantly  better  than  the  adaptive  controller.  In 
Fig.  8,  there  is  a  noticeable  overshoot  because  of  the  process 
gain  uncertainties,  whereas  the  MPC  controller  drops  the 
power  density  to  the  required  level  in  only  three  seconds,  de¬ 
spite  the  model  uncertainties.  It  is  important  to  note  that  the 
observed  phase  lag  in  the  MPC  response  in  due  to  the  fact 
that  the  solution  to  the  MPC  optimization  problem  is  only 
implemented  at  the  next  control  step.  Fig.  15  demonstrates 
the  effectiveness  of  MPC  when  operating  in  the  vicinity  of 
the  sign  change  in  static  gain.  Recall  that  the  sign  change 
for  the  modified  conditions  occurs  at  0.85  A/cm2.  As  seen 
in  Fig.  15,  the  control  variable,  the  current  density,  crosses 
that  value  with  no  ill  effects. 

Figs.  16  and  17  demonstrate  the  multivariable  capabilities 
of  MPC.  Here,  both  the  current  density  and  the  coolant  inlet 


time  [sec] 


Fig.  16.  Multivariable  control  using  the  current  density  and  coolant  inlet 
temperature. 


time  [sec] 


Fig.  17.  Multivariable  control  for  Pscl  =  0.6W/cm2. 


temperature  are  used  as  control  variables.  Fig.  16  shows  the 
performance  to  the  same  set  point  changes  as  seen  in  Fig.  15, 
in  which  only  current  density  is  used.  Comparing  the  results, 
it  can  be  concluded  that  the  addition  of  the  coolant  tem¬ 
perature  increases  the  complexity  of  the  controller  without 
contributing  to  the  performance.  On  the  other  hand.  Fig.  17 
shows  that  the  coolant  temperature  can  be  used  to  enable  the 
system  to  reach  values  that  would  otherwise  be  out  of  reach. 
Fig.  14  shows  that  the  current  density  alone  can  not  force  the 
system  to  reach  0.6  W/cm2,  whereas  using  both  variables,  as 
in  Fig.  17,  the  system  reaches  the  set  point  in  five  seconds. 


5.  Conclusions 

A  detailed  model  of  a  fuel  cell  and  its  simulation  has  iden¬ 
tified  severe  nonlinearities  in  the  response  to  a  change  in  the 
current  density.  The  sign  change  of  the  static  gain,  which 
occurs  within  the  normal  operating  range  of  the  device,  to¬ 
gether  with  the  uncertainty  of  the  precise  location  at  which 
this  occurs,  indicate  that  nonlinear  control  is  required  to  ad¬ 
equately  regulate  the  power  output  of  the  fuel  cell  in  the 
case  of  drastic  load  changes.  Use  of  nonlinear  model  pre¬ 
dictive  control  enables  accurate  control  in  the  case  of  such 
uncertainties,  with  multivariable  control  improving  perfor¬ 
mance.  The  use  of  the  MPC  algorithm  for  fuel  efficiency 
while  satisfying  load-change  demands  is  work  in  progress. 
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Appendix  A.  Derivation  of  the  electrochemical 
equations 


The  voltage  of  the  cell  as  a  function  of  the  current  density 
can  be  evaluated  as  follows 


Vceii  =  Voc  -  /?(•*)  - 


Kx)tm 

CTm(x)  ’ 


(39) 


where  the  first  term  on  the  right  is  the  open  circuit  potential. 
Yi  and  Nguyen  assumed  a  constant  open  circuit  potential, 
since  their  analysis  was  under  isothermal  conditions.  How¬ 
ever,  since  this  equation  will  be  used  under  transient  condi¬ 
tions,  the  Nernst  equation  is  used  to  account  for  temperature 
and  species  compositions 


Kc  =  K 


(40) 


The  second  term  is  the  activation  overpotential.  It  is  de¬ 
rived  from  the  Butler- Volmer  equation  [2]  for  the  case  where 
the  second  expression  (anodic)  can  be  neglected  (high  cur¬ 
rent  densities/overpotentials) 


/=/() 


"Co,,  six)  (- aFr](x)\ 

Ch2,sW  f(l-a)Frix\ 
CH2,b(x)eXPl  RT  ) 


(41) 


In  fuel  cells,  the  overpotential  is  mainly  at  the  cathode  so 
the  Tafel  equation  can  be  used 


tl(x)  = 


RTS 

- —In 

0.5  aF 


Kx)  \ 
/o  Pq,  (x) ) 


F  \l0Po2(x)J 

(42) 


where  7°  is  the  exchange  current  density. 

Yi  and  Nguyen  used  the  bulk  oxygen  partial  pressure  as 
the  partial  pressure  at  the  surface 


keen  =  v 


^in  (  I(X)  \  _  1{x)tm 

F  n\IoPo2(x))  <xm(x) 


(43) 


This  should  be  changed,  since  the  concentration  overpo¬ 
tential  is  not  taken  into  account.  Assuming  constant  diffu¬ 
sion  coefficient  in  the  gas,  D ,  and  a  diffusion  layer  thickness 
S,  we  can  say  that  at  steady  state  conditions  the  oxygen  ap¬ 
proaching  the  surface  is  equal  to  the  oxygen  reacting 


—D 


d 2_Po2 

dy2 


=  0  =>■ 


dPp2 

dy 


=  const 


(44) 


However,  bulk  and  surface  pressures  can  be  defined.  Inte¬ 
grating  yields: 


Po-,  (x,  S )  —  Po2  sW 

Po2(x,  y)  =  Po2Ax)  +  2  \  2  y  (45) 


and 


d  Pp2  (x,  y )  _  Pp2,b(x)  -  Pq2,s(x) 
dy  5 


(46) 


Thus,  comparing  the  oxygen  approaching  and  being  con¬ 
sumed  at  the  surface 


-D 


Po2(x,  8)  -  Pq2Ax) 


Kx) 
"  4  F 


and  isolating 


Po2,sW  =  Po2,b(x)  ~ 


In  the  same  fashion 


8I(  x) 
AFD 


Ph2,s  W  =  Ph2, b(x)  - 


Ph20,sW  =  PH20,bW  + 


8I(x) 
2FDHl 
8l(x) 


2  P  Du2o 

This  is  inserted  into  Eq.  (43)  giving 

/ 


RT 

Kell  =  V0°c  + 


[P|l2,bW  ( 2Fna2  ) 


[ 


Po2,b(x)  - 


SI(x) 


AFD, 


(47) 


(48) 


(49) 


o2 


-i  0.5 


Fu20,b(x)  +  (  2FO„*o  ) 


RT 

—In 

F 


Kx) 


I(x)t„ 


(9) 


It  is  clear  that  as  the  current  increases,  the  cell  voltage 
decreases.  The  voltage  on  the  electrodes  must  be  constant 
owing  to  the  high  conductivities  of  the  electrodes. 


Appendix  B.  Solution  of  current  density 


Eq.  (9)  is  a  relatively  difficult  equation  to  solve  numeri¬ 
cally  for  I(x),  since  it  does  not  behave  nicely  at  low  num¬ 
bers  and  negative  numbers  are  not  acceptable.  Neglecting  the 
concentration  overpotential  results  in  an  easier  expression 
to  solve.  Taking  an  exponent  of  both  sides  and  rearranging 
we  can  define 


m  = 


Kx) 


lo 

KhKu 

Ph2o 


exp 

,5/2 


Kell  -  KUC  + 


Kx)tm 

<7m(x) 


2  F 
RT 


=  0 


(50) 


Each  point  along  the  channel  must  satisfy  this  equation. 
This  is  a  relatively  easy  function  to  solve,  since  its  first  and 
second  derivatives  are  unconditionally  positive 
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d  /  .  1 

—  =2 — exp 

d  /  Io 


Kell  -  Kc  + 


I(x)tm  \ 
ffmW  / 


X 


Kx)  +  I(x)2 


tm  F  ' 
crm(x)  RT 


(51) 


This  type  of  function  is  well-suited  to  Newton’s  method, 
which  converges  very  quickly.  In  any  case  the  solver  backs 
up  Newton’s  method  with  bisection  iterations. 

If  the  concentration  overpotential  is  not  neglected,  Eq.  (9) 
is  manipulated  to  give 


f(J(x))  =  I  Pn2. bW 
I{x)2 

T- 
i0 


SI(x) 

~  2FDH2 

A^O.bto  + 
x  exp  (  (  Vcell  -  y0°c  + 


^Oo.bW  - 
SI(x)  \ 


SI(x ) 
4/-7J(>, 


2.5 


2FDu2oJ 

I(x)tm 


crm(x) 


)%) 


—  =0 


(53) 


whose  first  derivative 


d / 
d  7 


I(x)2S  exp(2 F/RT(Vcen  -  V°c  +  (7(x)rm/gm(x)))) 
2DHl0Fll 

2 1(x)  exp(2 F/RT(Vceil  -  V°c  +  (. I(x)tm/am(x )))) 
(4o,bW  +  (8I(x)/2FDH2o)) 


0.625 8(PH2,b(x)  ~  (8I(x)/2FDHJ)(P02,b(x) 
-(8I(x)/4FDo2))L5 


Dq2F 

8(P02,b(x)  -  (8I(x) / 4FDq2))2S 
2Dn2F 

2 1{x)2  exp(2F/RT(Vceu  -  V0°c  +  (Kx)tm/am(x)))) 
Ftm(Pn2o,b(x)  +  ( 8I(x)/2FDh2o )) 

Io2RTam(x ) 

(54) 


is  unconditionally  negative.  The  second  derivative,  however, 
may  change  signs,  which  can  be  problematic  for  Newton’s 
method  of  solution.  Consequently,  Newton’s  method  may 
experience  convergence  problems,  depending  on  the  values 
in  the  equation,  and  should  be  protected  using  bisection 
iterations. 

As  said  before,  this  value  of  7(x)  is  used  in  the  ODEs. 
Since  the  membrane  conductivity,  a,  temperatures  and  con¬ 
centrations  change  along  the  channel,  the  local  current  den¬ 
sity  will  vary,  as  well.  Note  that  since  all  the  changes  are 
continuous,  the  previous  solution  for  Kx)  is  an  excellent  ini¬ 
tial  guess  for  I(x  +  clx). 


Appendix  C.  Resolving  problems  when  integrating 
along  the  channel 

As  stated  in  Section  2,  numerical  solution  of  the  system 
equations  presents  several  challenges.  Firstly,  the  system 


equations  as  formulated  depend  on  whether  or  not  liquid 
water  is  present  in  the  channel.  Obviously,  if  liquid  water 
is  not  present  there  can  be  no  evaporation  regardless  of  the 
vapor  partial  pressure.  The  simple  solution  of  adding  a  logi¬ 
cal  switching  condition  introduces  stiffness  to  the  ODEs  be¬ 
ing  solved,  which  can  easily  lead  to  unstable  or  oscillatory 
numerical  solutions.  Instead,  the  ODE  solver  implemented 
uses  MATLAB’s  “events”  option,  which  checks  for  a  num¬ 
ber  of  situations: 


1.  Cathode  or  anode  liquid  reaches  zero. 

2.  Cathode  or  anode  vapor  partial  pressure  reaches  its  sat¬ 
urated  pressure. 


In  each  case,  the  ODE  solver  terminates  the  integration, 
returning  the  calculations  up  to  that  point  and  restarts  cal¬ 
culation  with  a  different  set  of  ODEs.  For  example,  if  there 
is  no  liquid  water  in  the  cathode  and  the  vapor  pressure  is 
below  saturation,  Eq.  (3)  is  not  used.  When  the  cathode  be¬ 
comes  saturated,  the  integration  terminates  and  is  restarted 
from  the  same  point,  this  time  using  Eq.  (3). 

In  addition,  it  is  noted  that  the  system  of  equations  in¬ 
volves  variables  with  scales  of  several  orders  of  magnitude 
apart-namely,  the  temperatures  are  in  K  of  o(102)  while  mo¬ 
lar  flow  rates  are  in  mol/s  of  o(10-6),  which  again  leads 
to  problems  in  their  numerical  solution.  This  is  resolved  by 
scaling  all  molar  flowrates  using  the  value  of  the  inlet  molar 
flow  of  hydrogen,  and  all  temperatures  using  the  tempera¬ 
ture  of  the  solid  at  x  =  0,  for  example 


A7h2(x)  =  Mh2(0)Mh2,  scaled  (x) 
T&{x)  =  Ts(0)T^  scaled  (A) 


The  ODEs  used  in  the  actual  integration  are  modified  to 
use  the  scaled  variables. 


Appendix  D.  Solution  of  transient  PDE 

The  transient  PDE  describing  the  temperature  profile  of 
the  solid  in  the  channel  direction  (Eq.  (18))  is  solved  using 
the  Crank-Nicholson  method.  This  method  approximates  the 
partial  differentials  by  finite  differences,  resulting  in  a  sys¬ 
tem  of  algebraic  equations,  which  describe  the  dependence 
of  the  temperature  profile  at  time  j  +  1  on  the  profile  at  time 
j.  This  system  is  solved  starting  at  the  initial  conditions,  with 
each  solution  giving  the  temperature  profile  at  the  next  time 
step. 

Introducing  the  following  expressions 


Y  =  P\  (7a  +  Tc)  +  fijTvj 

M  —  4 -  —  Mx 

m  —  iww,a,i+l  ™w,a,i-l  +  ^w,c,i+l  ™w,c,i-l 


(55) 
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Eq.  (17)  becomes 

7s,i,j+l  ~  Ts.i.j  _  7s,i+l,j+l/2  ~  2rs,i,j+1/2  +  rs,i—  1 0+1/2 


A  t 


Ax2 

(AH 

+  ^i.j+1/2  -  yT, s,i,j+i/2  -  e  ( 


VAcii 


X  7i,j+l/2  +  ^^vap(7i,j+l/2)M,j+l/2 


A  = 


case,  a  small  step  size  should  be  used,  to  ensure  stabil¬ 
ity.  Discretizing  Eq.  (19)  for  the  edges  and  rearranging 
gives 


—  (3k  +  2AxUc)TSiij+i  +  4kTsXj+i  —  kTs^j+i 
=  -2AxUcTinf 

kTsn—2,j+i  4kTsn—ij^.\  +  (3k  +  2AxUc)Tsnj^.\ 
=  2AxUcTinf 


(59) 


(56) 


Arranging  all  the  equations  into  matrix  form  and  defining 
matrices  A  and  B,  and  vector  C: 


—  (3k  +  2AxUc)  4  k  —k 


2A.y“ 

At 


0  •••  0 

0  •••  0 

+  2a  +  yAx2  —a  0  : 


(60) 


0 

0 


k  —4k  3k  +  2AxUc 


Inserting  expressions  for  the  solutions  on  the  “half-way 
point”  in  terms  of  values  on  the  solution  grid 


's,(+l. 7+1/2 


's,i,  7+1/2 


T’s.j-l,  7+1/2 


7s,  i+l,  7+1  +  7s,i+l,7 
2  ’ 
2s,  1,7+1  +  7s,  1,7 
2  ’ 

7s, /—1, 7+1  +  Tsj_\j 


(57) 


and  rearranging  gives 


B  = 


0  0  0 


0  a 


2A.y- 

At 


—  2a  —  yAx1 


0  0 

0  0 

a  0  : 


0  ■■  ■ 

0  •••  0 


a  a 

0  0  0 

(61) 


(2  Ax2  ,  . 

-aTSji+iJ+1  +  (  +2 a  +  yAx  )  TStiJ+l  -  aT!ij_L  j+1 


(  2  Ax2 

=  aTsj+ij  +  I  — - 2a  —  yAx')  TStiJ 

„  9  „  9  (AH 

+  aTSi^ij  +  2Ax  Yij+\/2  —  2Ax~e  I  +  Vceii 


X  Ii,j+l/2  +  ct>AxAHVilp(Tij+i/2)Mj  j+i/2 


(58) 


Recall  that  Y  and  I  are  calculations  based  on  Ts.  Solving 
this  entire  system  implicitly  requires  iterations  of  guessing 
the  profile  of  7's  at  time  j  +  1,  calculating  Y  and  I  at  time 
7+1,  and  inserting  them  into  the  equation.  The  right  hand 
side  will  then  be  a  vector,  and  the  linear  system  can  be  solved 
resulting  in  a  calculated  profile  for  Ts.  The  guess  for  7S  can 
then  be  iteratively  updated,  until  the  solution  for  Ts  at  time 
j  +  1  converges. 

Alternatively,  Y  and  7  can  be  used  explicitly,  in  which 
case  the  linear  system  can  be  solved  immediately.  In  this 


'0 

'0 

*2.7 

h,j 

C  =  2  Ax2 

Yn-lJ 

-  2Ax2s  ^  +  Vcell) 

In- 1,/ 

0 

0 

-2AxUcTiaf 

'0 

0 

AH^(T2j)M2j 

+  (p  Ax 

0 

A77vap(7’„_i,/)M„_i.7 

_  2AxUcTm( 

_0 

(62) 


The  linear  system  can  then  be  written  as 
ATj+ 1  =  BTj  +  C  (63) 

Note  that  matrices  A  and  B  remain  constant  for  all  times 
and  only  C  needs  to  be  updated  from  time  step  to  time  step. 
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If  the  system  is  to  be  integrated  implicitly,  Eq.  (58)  be¬ 
comes 


— «7s,i+l,;+l  + 


/  2  Ax2 
\  At 


+  2a  +  yAx  Tsjj+l  -  aTs^hJ+  j 


/  2  Ax2 

=  aTSyi+ij  +  I  — - 2a-  yAx2)  TSyiJ 

+  arSji_ ij  +  Ax2YiJ+1  +  A x2YUj  -  A x2s 

(AH  \  .  (AH  \ 

X  (  ^cel1  )  hj+ 1  —  A.V  £  (  -) p  *cell  )  h,j 

d>Ax  d>  Ax 

H - ^ —  A  Hmp  ( 7i_  j+ 1 )  Mj  y  j+ 1  H - — —  A  Hmp  (  7),  )  Mi  j 


(64) 


In  which  case 


"0 

"0 

*2.7 

0  i 

( AH  \ 

h-  .7 

*n-l,7 

—  Ax2£  j 

(2F  +  Vcel1) 

_0 

_0 

-2AxUcTinf 

"0 

0 

d>  Ax 
+  2 

A7/vap(72j)M2j 

0 

A77vap(r„_i,j)M„_  \  j 

_2AxUcTmf 

_° 

(65) 
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*2,7+1 


D{Tj+i)  =  Aa-2 


1,7+1 

0 


-  Ax2£ 


/  AH 
\2F 


Keen 


"0 

-2AxUcTmi 

*2,7+ 1 

0 

+ 

/i— 1,7+1 

0 

_0 

_2AxUcTmf 

0 


A  //yap  (  7),  j+ 1 )  M,;  j+ 1 


A  7/yap  ( 7),  7+ 1 )  Miy  j+ 1 
0 


(66) 


The  nonlinear  set  of  equations  is 

ATsJ+i  +  P(TsJ+  i)  =  +  C(7Vj)  (67) 

This  is  solved  using  a  relaxed  successive  substitution 
method  where 


Tsj+ 1  -  /(71s, J+l)  -  A-](BTSyJ  +  C(Tsj)  ~  D(TsJ+ 1)) 

(68) 

This  can  be  accelerated  or  damped  using  Wegstein’s 
method. 
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